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Submesoscale  Resolution 
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Abstract — In  this  paper,  a  variational  method  for  removing  po¬ 
sitioning  errors  (PEs)  from  drifter  trajectories  is  proposed.  The 
technique  is  based  on  the  assumption  of  statistical  independence 
of  the  PEs  and  drifter  accelerations.  The  method  provides  a  real¬ 
istic  approximation  to  the  probability  density  function  of  the  accel¬ 
erations  while  keeping  the  difference  between  the  filtered  and  ob¬ 
served  trajectories  within  the  error  bars  of  the  positioning  noise. 
Performance  of  the  method  is  demonstrated  in  application  to  real 
data  acquired  during  the  Grand  Lagrangian  Deployment  (GLAD) 
experiment  in  the  Northern  Gulf  of  Mexico  in  2012. 

Index  Terms — Computers  and  information  processing/data  pro¬ 
cessing,  mathematics/filtering  algorithms,  optimization,  smoothing 
methods. 


I.  Introduction 

STUDIES  of  the  world  ocean  with  Lagrangian  floats  have 
a  long  history  [24].  Drifters  were  mostly  employed  for  the 
exploration  of  the  large-scale  [25],  [19],  [7],  [26]  andmesoscale 
[9],  [27]  circulations  with  typical  intervals  between  the  drifter 
position  fixes  r(tn )  ranging  from  1  h  to  ten  days. 

In  recent  years,  there  has  been  a  growing  interest  in  the  anal¬ 
ysis  of  the  Lagrangian  coherent  structures  and  submesoscale 
processes  at  the  ocean  surface  (e.g.,  [3],  [23],  [11],  and  [4]),  fu¬ 
eled  by  the  necessity  to  increase  the  efficiency  of  oil  spill  man¬ 
agement,  search-and-rescue  missions,  and  ecological  control. 

The  Lagrangian  studies  at  submesoscale  resolution  require 
position  sampling  interval  St  as  small  as  several  minutes,  which 
amplifies  noise  in  the  estimates  of  drifter  accelerations 

a(tn)  =  [»•(£„_ i)  -  2 r(tn)  +  r(tn+i)]6t~ 2  (1) 

to  a  level  comparable  with  the  accelerations  experienced  by  the 
drifters.  Assuming  stationarity  and  statistical  independence  of 


Manuscript  received  February  19,  2014;  revised  August  07,  2014;  accepted 
August  21,  2014.  Date  of  publication  October  21,  2014;  date  of  current  version 
July  10,  2015.  This  work  was  supported  by  a  grant  from  The  Gulf  of  Mexico 
Research  Initiative  through  the  Consortium  for  Advanced  Research  on  Trans¬ 
port  of  Hydrocarbon  in  the  Environment  (CARTHE). 

Associate  Editor:  J.  F.  Lynch. 

M.  Yaremchuk  is  with  the  Naval  Research  Laboratory,  Stennis  Space  Center, 
MS  39556  USA. 

E.  F.  Coelho  is  with  the  Marine  Sciences  Department,  University  of  Southern 
Mississippi,  Stennis  Space  Center,  MS  39522  USA  (e-mail:  emanuel.coelho.ctr. 
po@nrlssc.navy.mil). 

Digital  Object  Identifier  10.1109/JOE.2014.2353472 


the  positioning  errors  (PEs),  (1)  yields  the  following  estimate 
of  the  PE-induced  acceleration  error: 


c tp  =  V6creSt  2 .  (2) 

For  a  sampling  interval  St  =  5  min  and  a  positioning  error  ae  = 
1.5  m,  the  acceleration  error  is  crp  =  4  x  10-5  m/s2,  a  value 
comparable  with  the  typical  Coriolis  acceleration  of  a  water 
parcel  moving  at  the  speed  of  2  m/s  in  midlatitudes.  On  the 
other  hand,  a  typical  acceleration  of  submesoscale  motions  in 
the  surface  layer  estimated  as  U2  / L  has  similar  magnitude  for 
a  characteristic  velocity  U  and  horizontal  L  scales  of  0.2  m/s 
and  1  km,  respectively  (e.g.,  [28]).  At  the  same  time,  subme¬ 
soscale  phenomena  are  characterized  by  significant  horizontal 
intermittency,  because  their  growth  and  development  are  associ¬ 
ated  with  mechanisms  (frontogenesis,  mixed-layer  instabilities) 
requiring  specific  environmental  conditions  (strong  subsurface 
stratification  and/or  horizontal  velocity  shear).  The  intermittent 
nature  of  submesoscale  turbulence  implies  that  the  acceleration 
probability  density  function  (pdf)  may  have  long  tails  and  a 
sharp  peak  in  the  region  of  small  (a  <  10“ 5  m/s2)  accelera¬ 
tions  corresponding  to  the  geostrophically  balanced  flows  with 
typical  velocities  of  0. 1-0.2  m/s. 

Accurate  estimates  of  the  Lagrangian  accelerations  in  the 
geophysical  flows  are  important  for  the  model  skill  improve¬ 
ment  because  they  provide  information  on  the  actual  forces 
acting  within  the  fluid.  In  addition,  the  Lagrangian  point  of  view 
is  in  many  aspects  the  most  natural  way  to  obtain  understanding 
of  the  turbulent  transport  and  mixing.  For  these  reasons,  there 
has  been  a  substantial  progress  in  both  experimental  [13],  [31], 
[18]  and  theoretical  [12],  [5],  [16]  studies  of  the  submesoscale 
structures  in  the  past  decade. 

Currently,  high-resolution  satellite  images  provide  the  best 
observational  window  into  the  submesoscale  world.  Devel¬ 
opments  of  the  towed  platforms  and  autonomous  underwater 
vehicle  technology  also  contribute  to  the  rapid  increase  of 
the  experimental  data  at  submesoscale.  In  the  recent  work  of 
Gaultier  et  al.  [10],  the  first  attempt  has  been  made  to  correct 
surface  velocity  field  using  submesoscale  sea  surface  tem¬ 
perature.  Drifter  trajectories,  sampled  at  temporal  resolution 
St  ~  L/U  <C  90  min,  constitute  another  valuable  source  of 
information  on  submesoscale  dynamics.  The  recent  Grand 
Lagrangian  Deployment  (GLAD)  experiment  conducted  in  the 
Gulf  of  Mexico  in  2012  appears  to  be  the  first  experimental 
study  of  submesoscale  processes  using  a  massive  array  of 
several  hundred  Lagrangian  drifters. 
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In  this  paper,  we  propose  a  variational  algorithm  for  filtering 
measurement  errors  from  drifter  coordinates  acquired  at  tem¬ 
poral  resolution  high  enough  to  contaminate  the  statistics  of 
accelerations  caused  by  physical  processes  in  the  upper  ocean. 
Since  the  typical  PE  magnitude  ae  is  small  compared  to  the  dis¬ 
tance  between  the  drifter  position  fixes,  PEs  have  a  small  impact 
on  the  Lagrangian  velocity  statistics,  but  distort  considerably  the 
statistics  of  Lagrangian  accelerations.  For  this  reason,  the  pro¬ 
posed  filtering  method  is  based  on  the  assumption  that  PEs  are 
statistically  independent  of  the  Lagrangian  accelerations  whose 
pdf  is  employed  for  regularization  of  the  cost  function. 

The  paper  is  organized  as  follows.  In  Sections  II  and  III, 
we  describe  the  method  and  the  data  used  for  its  validation.  In 
Section  IV,  the  method  is  applied  to  drifter  trajectories  docu¬ 
mented  during  the  GLAD  experiment.  The  results  are  summa¬ 
rized  and  discussed  in  Section  V. 


II.  Method 

Assuming  that  drifter  accelerations  a  and  errors  e  in  mea¬ 
suring  drifter  coordinates  are  statistically  independent,  their 
joint  (prior)  pdf  P(a,  e)  is 


P(a,e)  =  Pa(a)Pe(e).  (3) 

Consider  N  random  2d  vectors  an  and  en,  n  =  1, . . . ,  N, 
and  assume  that  Pe  can  be  approximated  by  the  Gaussian  dis¬ 
tribution,  so  that 


Pe{e  1,  •  •  •  ,eN)  oc  exp 


1 

2 


N 


E 


(4) 


where  crn  are  the  root  mean  square  (rms)  PE  variances.  In  ad¬ 
dition,  we  assume  that  the  prior  pdf  for  accelerations  Pa  can  be 
represented  by  the  stretched  exponential  law,  a  particular  case  of 
the  ^-exponential  distribution  with  q  =  1,  often  used  to  describe 
statistics  of  the  Lagrangian  accelerations  in  turbulent  flows  (e.g., 
[1]) 


Pa{*l,  •••,««)  OC  exp 


v 

hr  ^ 


(5) 


where  p  is  the  adjustable  parameter  defined  in  Section  IV. 

Let  us  also  assume  that  additional  information  on  a  and  e  is 
available  in  the  form  of  the  observed  drifter  positions  r*  and 
the  deterministic  relationships  (1)  and  en  =  rn  —  r* ,  with  rn 
being  the  unknown  true  position  at  time  tn.  Substitutions  of 
these  constraints  into  (3)  transforms  it  into  the  posterior  pdf  P* 
written  in  terms  of  r 


P*(r) 
oc  exp 


N- 1 
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(jn+l  -  2rn  +rK- i)2m 
(St2a)2^ 


N 
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2crn 


(6) 


The  pdf  (6)  can  be  also  obtained  using  the  well-known 
Bayesian  approach  under  the  standard  assumption  of  an  im¬ 
proper  prior  (e.g.,  [14]). 


The  most  likely  drifter  trajectory  can  be  found  by  maximizing 
the  probability  P*  [r]  in  the  2iV-dimensional  space  of  unknown 
true  positions  rn,  or  by  minimizing  the  following  cost  function: 

N—l  N 

J =  E  W° (r«+i  - 2rn  +r„_ ! )2^  +  Wn (r„  ■ ~K)2  -»■  min 

n= 2  n= 1 

(7) 

where  Wn  =  crn2/ 2  and  Wa  =  ( St2a )  2/C 

Since  the  true  acceleration  variance  cr2  is  unknown,  the  value 
of  Wa  is  still  undefined.  However,  it  can  be  expressed  in  terms  of 
the  directly  measured  quantities:  observed  acceleration  variance 
cr2  and  the  acceleration  variance  cr2  induced  by  PEs  (2).  Because 
PEs  and  accelerations  are  statistically  independent,  then 

=  er2  +  cr2  (8) 

and,  combining  (8)  with  the  linear  relationship  between  cr2  and 
a2,  one  can  obtain 

Wa  =  (St2a)-2»  =  [CSt4(a2  -  a2p)}~ »  (9) 

where  C  is  the  proportionality  coefficient  (see  the  Appendix). 

Given  the  statistical  model  (6),  the  minimization  problem  (7) 
and  (9)  has  the  only  free  parameter  fi ,  which  can  be  optimized  by 
choosing  its  value  such  that  the  rms  difference  between  the  most 
likely  and  observed  drifter  trajectory  is  equal  to  the  magnitude 
of  the  Global  Positioning  System  (GPS)  noise. 

Numerically,  the  cost  function  (7)  is  minimized  using  the 
double-precision  version  of  the  quasi-Newtonian  descent  algo¬ 
rithm  of  Gilbert  and  Lemarechal  [9].  The  high  machine  preci¬ 
sion  is  essential  because  trajectory  lengths  are  several  orders 
in  magnitude  larger  than  the  GPS  errors.  In  the  real  application 
described  below,  the  algorithm  converged  rather  quickly  and  re¬ 
quired  less  than  a  hundred  iterations  to  process  a  trajectory.  The 
total  central  processing  unit  (CPU)  time  needed  for  GLAD  data 
filtering  (few  million  position  fixes)  never  exceeded  2  min  on  a 
laptop. 


III.  Data 

A.  The  GLAD  Experiment 

In  summer/fall  2012,  an  extensive  Lagrangian  experiment 
was  conducted  in  the  northern  Gulf  of  Mexico.  An  array  of 
319  surface  drifters  was  deployed  in  late  July  2012,  and  the 
majority  of  drifters  persisted  throughout  October  2012,  a  lim¬ 
itation  imposed  by  their  battery  life.  The  deployed  SPAR-type 
drifters  were  different  from  those  traditionally  used  as  part  of 
the  global  drifter  program.  They  were  similar  to  the  surface 
drifters  (e.g.,  [8]  and  [22])  of  the  Coastal  Ocean  Dynamics  Ex¬ 
periment  (CODE),  originally  designed  to  follow  the  near  sur¬ 
face  upper-ocean  flows  in  the  presence  of  wind  and  waves.  As 
such,  the  GLAD  drifters  were  drogued  at  a  depth  of  approxi¬ 
mately  1  m  to  decouple  their  motion  from  direct  wind  forcing 
(estimated  l%-3%  windage  [9])  and  damp-wave-induced  mo¬ 
tions  without  introducing  a  wave-phase-related  bias  [9].  Posi¬ 
tions  were  reported  at  5 -min  intervals,  although  significant  data 
gaps  occurred  mostly  due  to  adverse  weather  conditions,  which 
required  careful  data  processing  to  produce  accurate  trajectory 
and  velocity  estimates. 
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Fig.  1.  Trajectories  of  the  GLAD  drifters  (shown  in  gray,  every  second  trajectory  is  given).  The  black  circle  depicts  the  location  of  the  Deep  Water  Horizon  oil 
spill.  The  distribution  of  the  time  intervals  between  the  drifter  position  fixes  is  given  in  the  inset. 


The  experiment  constituted  an  essential  first  step  to  study 
pollutant  transport  by  surface  currents  and  explore  interactions 
between  mesoscale  and  submesoscale  processes  (Fig.  1).  The 
major  objective  of  the  study  was  to  improve  understanding  of 
the  surface  transport  phenomena  in  the  open  ocean  and  near¬ 
coastal  regions  with  an  ultimate  goal  of  improving  Lagrangian 
predictability  and  dispersion  of  pollutants  at  time  scales  from  a 
few  hours  to  several  weeks. 

To  address  such  a  wide  spectrum  of  temporal  scales,  the 
drifters  were  equipped  with  Secure  POsiTioning  (SPOT, 
http://www.findmespot.com)  devices  that  recorded  drifter 
positions  every  St  =  316  s  on  average.  Overall,  more  than 
5  x  106  position  fixes  were  registered  with  the  total  duration  of 
the  drifter  trajectories  of  49.2  years.  The  drifters  were  drogued 
at  the  average  depth  of  1  m  and  were  released  in  a  variety  of 
clusters  to  explore  a  wide  range  of  spatial  scales  (see  [2 1  ]  and 
www.carthe.org/research.php  for  details),  tn  this  paper,  we 
focus  on  the  analysis  of  drifter  trajectories  with  an  objective  to 


remove  SPOT  positioning  errors  and  improve  the  estimates  of 
the  drifters’  accelerations. 

Although  the  time  intervals  between  the  position  fixes  were 
not  constant,  the  majority  of  them  (85%)  were  within  several 
seconds  of  St  =  5  min  (inset  in  Fig.  1),  with  the  multiples  of 
St  corresponding  to  another  14%  of  the  total  number.  The  av¬ 
erage  distance  between  the  sequential  position  fixes  was  close 
to  1 1 5  m,  providing  the  mean  drifter  velocity  of  0.36  m/s.  These 
numbers  show  that  even  at  such  a  short  sampling  time  interval, 
PEs  of  several  meters  introduce  velocity  errors  of  a  few  centime¬ 
ters  per  second,  which  is  acceptable  for  most  of  the  applications, 
tn  contrast,  errors  in  drifter  accelerations  appear  to  be  quite  sub¬ 
stantial  [see  (2)]  with  the  magnitude  several  times  larger  than  the 
typical  magnitude  of  water  acceleration  in  surface  layer.  Despite 
the  high  level  of  acceleration  noise  introduced  by  the  GPS  er¬ 
rors  at  such  a  high  temporal  resolution,  the  task  of  filtering  this 
noise  becomes  feasible  because  the  PE  statistics  can  be  assessed 
with  a  reasonable  degree  of  accuracy. 
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To  account  for  these  tails,  we  elected  a  nonlinear  approach: 
if  the  trajectory  r(t)  experienced  a  sudden  excursion  causing  a 
large  spike  in  acceleration,  the  PE  error  variance  was  inflated  as 
follows: 

2 

an  =  <7e[l  +ll9{ln  ~  1)],  7 'n=~J  (10) 

ac 

where  ac  is  the  cutoff  acceleration  and  0  is  the  step  function. 
In  actual  computations,  the  relationship  between  a  and  r  was 
modified  to  account  for  variations  of  the  time  intervals  between 
the  position  fixes 

_  (rn+i  —  rn)^tn1  —  (rn  —  rn-i)5tn_1  /V1  ^ 

an  —  (r,  I  e  I  \  •  (II) 

yOtn  T  otn— i  j 
2 


Fig.  2.  Scatter  plot  of  the  coordinates  of  the  drifter  309  while  it  was  still  on 
shore.  Approximation  of  its  meridional  PE  by  the  Gaussian  pdf  is  shown  in  the 
inset. 


Fig.  3.  Relationships  between  cra  and  ac  used  for  identification  of  the  cutoff 
acceleration.  The  acceleration  values  (ms-2)  are  multiplied  by  105. 


B.  Positioning  Errors 

Among  more  than  five  million  position  fixes  available  in  the 
GLAD  database,  a  certain  number  were  made  from  drifters  that 
were  not  yet  deployed  did  not  change  their  location.  In  partic¬ 
ular,  the  stationary  position  of  the  drifter  309  was  recorded  for 
more  than  two  weeks  (4980  position  fixes)  while  it  was  on  board 
of  the  docked  deployment  vessel  at  29.3886  °N  90.7145  °W. 
This  allowed  us  to  assess  the  statistics  of  PEs  under  the  undis¬ 
turbed  ashore  conditions  (Fig.  2). 

In  the  first  approximation,  the  PE  distribution  can  be  de¬ 
scribed  by  the  uncorrelated  Gaussian  pdfs  with  the  variances 
of  a x  =2.1  and  ay  =  1.9  m  in  the  zonal  and  meridional 
directions,  respectively  ( ae  =  y/cr^  +  a^  =  2.8  m).  The  PE-in- 
duced  “acceleration”  rms  variance  ap  was  computed  using  the 
recorded  time  variation  of  the  drifter's  coordinates  and  was 
found  to  be  7.6  x  10-5  ms-2. 

It  is  noticeable  that  the  pdf  in  Fig.  2  exhibits  significant  non- 
Gaussian  tails  for  the  errors  exceeding  5  m.  Accurate  fitting  of 
these  tails  by  commonly  used  analytic  functions  (e.g.,  stretched 
exponentials)  is  impossible  because  of  the  insufficient  statistics. 
Moreover,  these  tails  are  likely  caused  by  random  events  (e.g., 
a  sudden  persistent  decrease  in  the  number  of  satellites)  whose 
statistics  are  difficult  to  model. 


The  parameter  ac  can  be  computed  by  taking  into  account 
the  statistical  independence  of  the  PEs  and  drifter  accelerations 
and  considering  limitations  on  the  estimates  of  the  acceleration 
probability  density  in  the  tails,  imposed  by  the  finite  volume  of 
the  GLAD  data.  First,  aa  can  be  estimated  from  (8),  (10),  and 
(11)  as  a  function  of  ac,  given  the  observed  values  of  cr*  = 
8.9  x  10-5  ms-2  and  ap  =  7.6  x  10-5  ms-2.  With  the  de¬ 
crease  of  ac,  the  value  of  cr2  increases  as  more  points  on  the 
trajectory  are  treated  as  outliers  and  penalized  by  large  obser¬ 
vation  variance  according  to  (10).  As  a  consequence,  the  value 
of  a a  =  ^ al  —  crp(ac )  decreases  until  it  reaches  zero  at  ac  = 

1.39  x  10-4ms-2.  This  dependence  is  shown  by  the  solid  line 
in  Fig.  3.  The  second  relationship  between  aa  and  ac  can  be 
obtained  by  considering  the  limitation  in  the  assessment  of  the 
acceleration  statistics:  With  the  finite  number  of  samples  M,  the 
remote  (a  ac)  structure  of  the  pdf  tails  cannot  be  adequately 
resolved  if  these  tails  are  described  by  less  than  a  few  hundred 
samples,  i.e.,  MPa(a  >  ac)  ^  102.  Using  the  analytic  pdfs  de¬ 
scribing  the  experimental  data  on  the  Lagrangian  accelerations, 
e.g.,  [1]  and  [2],  it  is  possible  to  estimate  that  for  several  million 
GLAD  samples  the  cutoff  value  can  be  reasonably  constrained 
by  ac  =  6a a  (a  dashed  line  in  Fig.  3).  Among  the  two  roots 
visible  in  Fig.  3,  we  used  the  smaller  one,  primarily  because  the 
larger  cutoff  value  identified  too  few  spikes  as  outliers. 

One  should  note  that  for  a  deployed  drifter  the  value  of  ap 
should  be  larger  than  7.6  x  10-5  m/s2  because  of  the  addi¬ 
tional  factors  causing  the  loss  of  sight  of  GPS  satellites.  How¬ 
ever,  in  view  of  (8),  the  (hardly  measurable)  increase  in  ap  is 
constrained  from  above  by  the  value  of  the  observed  drifter  ac¬ 
celeration  a * .  In  the  case  of  GLAD  data,  this  increase  can  be 
considered  as  negligible  in  the  first  approximation  because  the 
value  of  aa  =  (cr2  -  cr2)1/2  =  4.5  x  10-5  m/s2  appears  to  be 
fairly  consistent  with  the  typical  accelerations  of  the  drifters  in¬ 
volved  in  submesoscale  motions. 

As  a  result  of  PE  analysis  of  the  GLAD  data,  the  cost  function 
was  modified  to  account  for  non-Gaussian  tails  in  the  PE  distri¬ 
bution.  The  modification  allows  PE  variances  to  change  in  time 
through  parameterization  (10)  with  the  cutoff  scale  ac  =  1 .56  x 
10-4  ms-2,  corresponding  to  the  Coriolis  acceleration  experi¬ 
enced  by  a  water  parcel  traveling  at  a  speed  of  2.2  m/s.  This 
value  corresponds  to  the  Rossby  number  Ro  =  U/  fL  ~  30, 
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Fig.  4.  Time  evolution  of  the  observed  (black)  and  filtered  (light  blue)  acceleration  magnitude  of  a  GLAD  drifter  trajectory.  Yellow  rectangle  shows  time  interval 
corresponding  to  the  piece  of  the  trajectory  exposed  in  Fig.  5(b)-(c).  The  magnitude  of  accelerations  is  divided  by  the  mean  Corioilis  parameter  /  =  7.3  x 
IQ"5  s"1. 


which  is  far  beyond  the  typical  range  0.5  <  Ro  <  5  of  subme- 
soscale  motions  and,  therefore,  has  a  negligible  impact  on  the 
acceleration  statistics  of  submesoscale  turbulence. 

IV.  Results 

Before  applying  the  filter  to  GLAD  trajectories,  the  raw  data 
were  quality  controlled.  First,  12-h  pieces  of  the  trajectories 
were  cut  from  both  ends  to  remove  time  intervals  when  the 
drifters  were  on  board  before  and  immediately  after  the  deploy¬ 
ment.  The  trajectories  with  the  maximum  speeds  exceeding  3 
m/s  and  shorter  than  five  days  were  also  discarded.  Finally,  the 
trajectories  with  the  maximum  sampling  time  interval  larger 
than  2  h  and  persistent  (lasting  more  than  2 St)  zero  accelera¬ 
tions  were  discarded  to  remove  the  impact  of  time  extrapolation 
of  the  SPOT  software  in  case  when  satellite  signal  was  lost. 

The  filter  was  tuned  by  performing  a  series  of  filtering  exper¬ 
iments  with  variable  parameter  fi.  Statistically,  the  mean  rms 
difference  ae  =  2.8  m  between  the  raw  and  filtered  drifter  coor¬ 
dinates  was  achieved  at  /i  =  0.91  and  this  value  was  chosen  as 
an  optimal  one  for  the  filter.  Interestingly,  it  is  close  to  the  one 
(0.90)  obtained  by  Lewis  and  Swinney  [15],  for  the  distribution 
of  the  Lagrangian  velocity  increments  in  the  turbulent  flow. 

In  total,  169  trajectories  with  2.4  million  position  fixes  were 
filtered.  Among  those,  51  387  acceleration  spikes  were  identi¬ 
fied  and  92%  of  them  were  reduced  to  a  level  below  ac  after  the 
PE  correction  of  the  trajectories.  The  remaining  accelerations 
never  exceeded  1.3ac. 

A  comparison  of  the  raw  and  filtered  trajectories  is  shown  in 
terms  of  accelerations  (Fig.  4)  and  actual  coordinates  (Fig.  5)  for 
an  arbitrarily  selected  drifter.  During  2.5  months  of  the  drifter 
operation,  there  were  multiple  events  of  enhanced  accelerations 
lasting  from  several  minutes  to  several  days  (Fig.  4).  Of  course, 
the  origin  of  accelerations  cannot  solely  be  attributed  to  the  GPS 
noise,  because  the  drifters  are  also  exposed  to  a  wide  range  of 
unresolved  small-scale  phenomena  such  as  sudden  wind  gusts, 
Langmuir  cells,  wave-induced  Stokes  drift,  etc.  However,  ac¬ 
cording  to  the  estimates  in  Section  III,  more  than  half  of  the 
acceleration  variance  (visible  in  black  in  Fig.  4)  should  be  at¬ 
tributed  to  the  GPS  noise. 

Fig.  5  provides  a  more  detailed  exposure  of  the  filter's  cor¬ 
rection  of  the  drifter's  coordinates.  Of  particular  interest  are  the 


abrupt  changes  between  the  periods  of  relatively  uniform  mo¬ 
tion  and  periods  of  intense  variations  of  the  drifter  position  at 
the  scales  of  less  than  a  few  kilometers  [Fig.  5(a)].  Such  be¬ 
havior  is  indicative  of  intermittent  nature  of  the  submesoscale 
turbulence.  The  part  of  the  trajectory,  corresponding  to  the  se¬ 
quence  of  high  acceleration  events  on  day  38  visible  in  Fig.  4,  is 
shown  in  Fig.  5(b).  Because  the  resolution  of  Fig.  5(b)  is  still  in¬ 
sufficient,  Fig.  5(c)  shows  deviations  of  the  trajectory  from  the 
uniform  motion  between  its  end  points.  An  example  of  filtering 
the  sudden  coordinate  excursion  observed  on  day  42  (see  accel¬ 
eration  peak  exceeding  10/  m/s2)  is  demonstrated  in  Fig.  5(d). 
This  event  is  likely  to  be  caused  by  the  GPS  error  because  the 
drifter  immediately  returned  to  the  original  track.  Random  in¬ 
spection  of  several  events  of  this  type  has  shown  that  the  filter 
performed  reasonably  well. 

The  blue  curve  in  Fig.  6  shows  the  pdf  of  the  accelerations 
of  the  analyzed  GLAD  trajectories.  In  the  region  \a\  <  0.9 ac, 
the  pdf  is  virtually  indistinguishable  from  the  Gaussian  noise 
parabola  with  the  rms  acceleration  error  variance  corresponding 
to  the  positioning  error  of  2  m.  This  supports  the  validity  of  the 
PE  model  used  in  the  formulation  of  the  cost  function  (7).  The 
clearly  visible  tails  of  the  observed  acceleration  pdf  stretch  to 
the  values  of  3.5  x  10-4  m/s2  and  could  be  attributed  to  the  non- 
Gaussianity  of  the  error  distribution  in  Fig.  2.  They  are  removed 
by  the  filtering  algorithm.  The  deviations  of  the  filtered  drifter 
positions  from  the  observed  ones  rarely  exceed  several  meters, 
and  appear  to  be  negligible  compared  to  the  spatial  scales  of  the 
submesoscale  motions  [Fig.  5(a)].  However,  the  pdf  of  filtered 
accelerations  has  changed  dramatically,  and  appears  to  be  more 
realistic. 

To  support  this  statement,  we  calculated  the  pdf  of  the 
Lagrangian  accelerations  from  the  1-km  Navy  Coastal  Ocean 
Model  (NCOM)  simulation  during  the  GLAD  experiment.  The 
virtual  trajectories  were  computed  by  integrating  the  NCOM 
velocity  field  contaminated  by  a  random  walk  process  whose 
effective  dispersion  coefficient  (150  m2/s)  was  specified  as 
the  mean  value  of  the  Smagorinsky  diffusion  coefficient  in 
the  model  run.  The  time  integration  was  performed  with  the 
fourth-order  Runge-Kutta  scheme.  Five  million  virtual  drifter 
positions  in  more  than  a  thousand  randomly  selected  trajecto¬ 
ries  were  sampled  at  5 -min  resolution  and  processed  to  obtain 
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dx,  m  x,  km 

(c)  (d) 

Fig.  5.  Trajectory  of  a  typical  GLAD  drifter  relative  to  its  mean  position  (a).  The  yellow  circle  shows  the  location  of  the  trajectory  part  expanded  in  (b).  Deviations 
of  the  trajectory  in  (b)  from  the  uniform  motion  between  the  end  points  are  shown  in  (c).  Yellow  circles  show  the  observed  positions,  and  the  blue  line  is  the  filtered 
trajectory.  Panel  (d)  is  the  same  as  in  (b),  but  for  the  high-acceleration  event  on  day  42. 


the  pdf  shown  in  Fig.  6  by  the  light  gray  line,  which  fits  the 
filtered  pdf  (the  red  line  in  Fig.  6)  in  a  reasonably  wide  range 
of  accelerations,  ft  is  also  remarkable  that  the  shape  of  the 
filtered  pdf  is  qualitatively  consistent  with  the  typical  pdfs  of 
the  Lagrangian  accelerations  observed  in  the  laboratory  exper¬ 
iments  (e.g.,  [17]  and  [30]  )  and  appears  to  be  in  a  reasonable 
quantitative  agreement  with  the  log-normal  superstatistics  of 
the  Lagrangian  accelerations  [2]  shown  by  the  dashed  line  in 
Fig.  6. 

The  impact  of  the  filter  on  the  correlation  functions  and 
spectra  of  the  Lagrangian  accelerations  is  shown  in  Fig.  7.  The 
curves  were  obtained  by  averaging  the  trace  of  the  autocorrela¬ 
tion  matrix  over  the  ensemble  of  the  drifter  trajectories  for  the 
delay  times  r  between  5  min  and  one  month.  The  most  severe 
contamination  of  the  acceleration  statistics  by  the  GPS  noise 
occurs  at  time  scales  close  to  the  sampling  interval,  rendering 
the  shape  of  sample  autocorrelation  to  be  close  to  the  5-function 
(a  light  gray  curve  in  the  inset  of  Fig.  7).  After  filtering,  the 
autocorrelation  exhibits  a  more  realistic  behavior,  better  visible 
in  spectral  representation:  In  contrast  to  the  raw  acceleration 
spectrum,  which  is  close  to  white  for  the  periods  T  <  5  h,  the 


spectrum  of  filtered  accelerations  exhibits  a  deep  minimum  at 
T  ~  10-15  min,  consistent  with  the  structure  of  wind  spectra 
(mesoscale  valley)  at  the  sea  surface  (e.g.,  [20]).  There  is  also 
a  slight  indication  of  an  upward  trend  at  T  <  10  min,  corre¬ 
sponding  to  the  atmospheric  microscale  maximum  observed  in 
the  3  <  T  <  0.5  min  range  (e.g.,  [29]).  The  high-frequency 
part  of  the  velocity  spectrum  is  affected  by  the  GPS  noise  to 
a  much  lesser  extent,  and  exhibits  a  structure  similar  to  the 
filtered  acceleration  spectrum  in  the  high-frequency  range. 

To  better  demonstrate  the  filter  performance,  we  processed 
the  GLAD  data  using  a  much  simpler  technique,  which  removes 
spikes  in  accelerations  by  prescribing  homogeneous  motion  on 
the  respective  time  intervals:  when  \an\  >  ac,  rn  in  (1 1)  was  re¬ 
placed  by  the  value  setting  the  numerator  to  zero.  Compared  to 
the  proposed  filter,  this  procedure  required  ten  times  less  com¬ 
puter  time  (12  s  for  processing  the  GLAD  trajectories),  but  pro¬ 
vided  moderate  improvement  of  the  acceleration  statistics  (dark 
gray  curves  in  Fig.  7):  the  autocorrelation  function  is  still  close 
to  the  5  shape,  while  the  spectrum  appears  to  be  homogeneously 
damped  at  high  frequencies.  The  resulting  pdf  (not  shown)  also 
demonstrates  an  unrealistic  behavior:  it  is  abruptly  cut  to  zero 
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Fig.  6.  Pdf  of  the  observed  GLAD  drifter  accelerations  (blue),  the  Gaussian 
model  of  the  accelerations  induced  by  the  GPS  noise  (black  parabola),  the  pdf 
of  the  filtered  accelerations  (red),  and  the  pdf  of  the  log-normal  superstatistics 
[2]  approximating  the  Lagrangian  accelerations  observed  in  the  laboratory  ex¬ 
periment  (the  dashed  line).  The  pdf  of  Lagrangian  accelerations  simulated  by  the 
NCOM  run  at  1-km  resolution  is  given  in  light  gray.  Thin  vertical  lines  show 
the  value  of  the  cutoff  acceleration. 


T,  min 

Fig.  7.  Spectra  of  the  observed  (gray)  and  filtered  (black)  Lagrangian  acceler¬ 
ations  averaged  over  the  ensemble  of  the  GLAD  trajectories.  Accelerations  are 
normalized  by  the  local  Coriolis  parameter  /  =  7.3  x  10-5  s_1 .  The  structure 
of  the  corresponding  autocorrelation  functions  near  the  origin  is  shown  in  the 
inset. 

at  \a\  >  ac  and  has  a  bulge  at  a  =  0,  which  is  the  result  of 
the  accumulation  of  the  filtered  spikes.  Of  course,  this  simple 
algorithm  can  be  improved  (e.g.,  by  introducing  spline  instead 
of  linear  interpolation),  but  these  changes  will  make  it  noncom¬ 
petitive  with  the  proposed  filter.  Besides,  we  believe  that  the 
computational  expense  of  the  proposed  filter  is  quite  modest  in 
atypical  application  involving  103  drifters  operating  during  102 
days. 

Statistical  consistency  with  the  NCOM  simulations  and  inde¬ 
pendent  experimental  data  indicates  that  the  filtered  trajectories 


are  substantially  more  realistic  in  terms  of  the  acceleration  sta¬ 
tistics.  At  the  same  time,  the  proposed  filter  has  little  impact  on 
the  velocity  statistics  because  deviations  of  the  filtered  trajecto¬ 
ries  from  the  observed  ones  are  relatively  small  and  introduce 
only  minor  (~0.5 — 1  cm/s)  corrections  to  the  drifter  velocities 
derived  from  the  differences  in  their  consecutive  positions. 

V.  Conclusion 

A  method  for  filtering  PEs  from  drifter  trajectories  is  pro¬ 
posed.  Its  performance  is  demonstrated  applying  it  to  the  GLAD 
data.  The  method  is  based  on  the  assumption  of  the  statistical 
independence  of  the  PEs  and  Lagrangian  accelerations.  The 
PE  statistics  are  approximated  by  the  ^-correlated  Gaussian 
process  with  time-dependent  variance  which  prescribes  higher 
PEs  to  outliers  and  identifies  them  by  spikes  in  the  acceleration 
time  series.  Statistics  of  drifter  accelerations  are  specified  by 
the  stretched  exponential  law,  whose  adjustable  parameter 
p  =  0.91  is  defined  by  optimizing  a  posteriori  PE  variance. 

Statistics  of  the  accelerations  obtained  after  filtering  the 
GLAD  data  are  consistent  with  the  results  of  high-resolution 
numerical  simulations  and  is  in  good  qualitative  agreement 
with  the  typical  pdf  shape  of  the  Lagrangian  accelerations 
observed  in  the  laboratory  experiments.  Quantitatively,  the 
obtained  pdf  is  similar  to  the  log-normal  superstatistics  of  Beck 
[2]  describing  the  distribution  of  the  Lagrangian  accelerations 
observed  in  a  developed  turbulent  flow  [32],  [33].  These  results 
indicate  potential  validity  of  the  filter  when  applied  to  other 
drifter  data  sets  acquired  at  high  temporal  resolution. 

It  is  noteworthy  that  a  simple  Gaussian  assumption  fi  —  1 
on  the  acceleration  statistics  produces  a  posteriori  pdf  with  a 
stronger  cutoff  of  the  tails  resembling  the  one  obtained  from 
the  NCOM  model  (the  gray  line  in  Fig.  6).  At  the  same  time,  a 
posteriori  rms  deviation  of  the  filtered  trajectories  from  the  ob¬ 
served  ones  was  significantly  different  from  the  observed  value 
ae  =  2.8  m,  because  the  Gaussian  distribution  appeared  to  be 
underdispersive.  Introduction  of  the  extra  tuning  parameter  p 
improved  the  situation  and  appears  to  be  reasonable  in  the  light 
of  the  cited  experimental  data  on  the  Lagrangian  accelerations. 

The  analysis  also  revealed  that,  on  average,  drifter  posi¬ 
tioning  errors  did  not  experience  significant  increase  from  their 
nominal  values  diagnosed  ashore.  This  can  be  partly  attributed 
to  the  relatively  calm  weather  during  the  GLAD  experiment  and 
the  sufficiently  long  GPS  signal  accumulation  (5  min)  between 
the  position  fixes  as  compared  to  the  time  scales  (5-20  s)  of  the 
wave-induced  disturbances. 

The  proposed  filtering  methodology  can  be  improved  in  sev¬ 
eral  ways.  First,  we  did  not  take  into  account  temporal  corre¬ 
lations  of  the  drifter  accelerations  along  the  trajectory.  Second, 
more  sophisticated  distributions,  such  as  ^-exponential,  or  log¬ 
normal  superstatistics  could  be  used  as  the  background  pdfs  for 
regularization  of  the  cost  function.  These  distributions  are  char¬ 
acterized  by  a  sharp  discontinuity  at  the  origin  ( a  =  0)  and  may 
be  difficult  to  handle  by  the  descent  algorithms  based  on  the  as¬ 
sumption  of  differentiability  of  the  cost  function.  In  fact,  during 
the  search  for  an  optimal  p,  we  did  encounter  a  significant  slow¬ 
down  of  the  convergence  rates  at  p  <  0.7.  Finally,  the  PE  distri¬ 
bution  can  also  be  improved  by  updating  the  GPS  noise  model: 
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the  PEs  are  likely  to  be  correlated,  at  least  during  the  sporadic 
events  of  signal  loss  from  a  significant  number  of  satellites. 

An  alternative  way  of  improving  the  quality  of  drifter  data  at 
high  temporal  resolution  is  to  reduce  the  PEs  themselves.  This 
can  be  easily  done  in  the  near-coastal  regions  by  using  the  ex¬ 
isting  network  of  the  continuously  operating  reference  stations 
(CORSs),  which  can  easily  deliver  positions  accurate  to  several 
centimeters.  In  the  open  ocean,  however,  this  method  may  not 
work  without  involving  a  sort  of  bottom-referenced  navigation 
and/or  more  sophisticated  real-time  ionosphere  correction  algo¬ 
rithms. 

We  consider  this  study  to  be  an  attempt  at  improving  the 
quality  of  the  drifter  data  acquired  at  high  temporal  resolution. 
The  methodology  can  be  extended  for  interpolation  of  the  drifter 
positions  on  a  regular  time  grid  by  introducing  the  local  time-in¬ 
terpolation  operators  in  the  first  (observational)  term  of  the  cost 
function.  There  is  also  plenty  of  room  for  further  development 
of  the  filter  along  the  aforementioned  lines. 


Appendix 

The  acceleration  variance  a2  is  given  by  the  integral 
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Introducing  polar  coordinates  ar  =  \a\  =  (a2  +  ajj)1/2  and 
Lp  =  tan ~1(ay/ax)  and  a  new  variable  rj  =  ( ar/a)2fl ,  the 
normalization  integral  is  reduced  to 


a 

2  (i 


p2n  p  oo 

/  d(p  r]1^~1e~r]dri 

Jo  Jo 


™'-r  (T)  =  »-■ 


(Al) 


where  the  gamma  function  is  defined  by  T(s)  = 
J0°°  xs~1e~xdx.  In  the  similar  manner,  the  first  integral  can 
be  rewritten  as 


r ^  f0°v2/tl-1e-ridfi  =  N— r(-V 

(A2) 

Substitution  of  N  from  the  right-hand  side  of  (Al)  into  the 
right-hand  side  of  (A2)  yields  the  expression  for  the  proportion¬ 
ality  coefficient  C  in  (9) 
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